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<N ; Abstract 
> . 

\ We simulate lattice QCD with 2 flavours of massless quarks on lattices of temporal extent A^^ = 8, 

■ to study the finite temperature transition from hadronic matter to a quark-gluon plasma. A mod- 

o ■ 

\^ . ified action which incorporates an irrelevant chiral 4-fermion interaction is used, which allows 



o 



simulations at zero quark mass. We obtain excellent fits of the chiral condensates to the magneti- 
zations of a 3-dimensional 0(2) spin model on lattices small enough to model the finite size effects. 

^ ■ This gives predictions for correlation lengths and chiral susceptibilities from the corresponding 

^ ■ 

spin-model quantities. These are in good agreement with our measurements over the relevant 

•i-H . 

^ ' range of parameters. Binder cumulants are measured, but the errors are too large to draw definite 

conclusions. From the properties of the 0(2) spin model on the relatively small lattices with which 
we fit our 'data', we can see why earlier attempts to fit staggered lattice data to leading-order 
infinite- volume scaling functions, as well as finite size scaling studies, failed and led to erroneous 
conclusions. 
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I. INTRODUCTION 



When hadronic matter is heated to temperatures of 150-200 MeV it makes a transition 
to a quark-gluon plasma. Such temperatures were present in the early universe. There is 
also evidence that temperatures this high can be observed in relativistic heavy-ion colliders 
such as RHIC. Indeed, claims that the quark-gluon plasma has been observed have been 
made at RHIC and at CERN. 

Lattice QCD enables one to simulate QCD at finite temperature and to examine the 
nature of this transition. When hadronic matter passes through this transition to the quark- 
gluon plasma, quarks and gluons are no longer confined. In addition, the approximate chiral 
symmetry associated with light quarks, which is broken spontaneously in cold hadronic 
matter, is restored. If we reduce the quark mass to zero, this transition must necessarily 
become a phase transition. We will consider the case of 2 light quarks {u and d), and take 
their masses to zero, assuming that the strange quark, which we neglect, is too heavy to 
induce a first-order transition. For 2 light quark fiavours, it is believed that the zero-mass 
transition is a second-order phase transition, which becomes a crossover with no true phase 
transition as soon as the quarks gain mass Q]. It is argued that the universality class of 
this phase transition is that of a 3-dimensional 0(4) spin model (Since the symmetry of 
massless 2-fiavour QCD is SU{2) x SU{2), which is isomorphic to 0(4), the tt-ct multiplet 
which provides the Goldstone pions, expresses 0(4) symmetry.) The reduced symmetry of 
staggered quarks reduces this symmetry to 0(2), so we might expect that lattice QCD with 2 
staggered quark fiavours has a phase transition in the universality class of the 3-dimensional 
0(2) spin model. 

In the standard formulations of lattice QCD, the Dirac operator becomes singular at zero 
quark mass, and so zero-mass simulations are impossible and very light quark simulations 
become very expensive. We have introduced a modified version of the staggered quark 
action, which has an additional irrelevant chiral 4-fermion interaction, "xQCD" Q]. Being 
irrelevant, this 4-fermion interaction does not affect the continuum limit, but it does render 
the Dirac operator non-singular, even in the chiral limit. Hence we simulate 2-flavour QCD 
with massless quarks using the xQCD action. 

We present the results of simulations of massless two-fiavour xQCD on Nt = 8 (12^ x 8, 
16^ X 8 and 24'^ x 8) lattices, in the neighbourhood of the finite-temperature transition. 
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State-of-the-art simulations with conventional and highly improved staggered-quark actions 
have been run at masses which are too large to accurately determine the universality class 
of the transition 0, Q. As a consequence, attempts to determine the universality class of 

;hat they are those for 



this transition, or to fit its scaling properties under the assumption t 
0(4) or 0(2) have produced confusing results Q, ^,1? J|liJlO, 11, 12|. 

Our earlier work on A^^^ = 4 2, and A''^ = 6 14. Iisl] lattices had produced results 
which appeared to preclude either 0(2) or 0(4) universality. We had hoped that running at 
Nf = 8 would broaden the transition sufficiently to allow a better determination of its nature. 
However, preliminary analysis of our measurements indicated that the chiral condensate 
decreases too fast for an 0(2) or 0(4) interpretation. The major problem appears to be 
that our 'data' needs to be extrapolated to infinite volume, and all such extrapolations 
became unreliable close to the transition. For this reason we abandoned our attempts to 
extrapolate to infinite volume, and instead fit our chiral condensates to the magnetization 
curves for the 3-dimensional 0(2) spin model, also on finite volumes. By choosing the lattice 
volumes for best fits, we are able to find good fits to our chiral condensates over most of the 
range of our measurements. We are also able to use our spin model fits to make acceptable 
predictions for the correlation lengths and chiral susceptibilities. 

Once we establish that the behaviour of lattice QCD appears to be consistent with that 
of the 0(2) spin model, we are able to use the properties of this simple model obtained both 
from our own simulations and the extensive literature on this model 0,0,0,0 
to understand why naive attempts to establish 0(2) universality failed. The sizes of the 
0(2) lattices that give good fits to our xQCD 'data' are too small (5^-9^) to be anywhere 
close to the infinite volume limit. In addition, the region over which the leading term suffices 
to describe the critical scalin g of the magnetization of the 0(2) spin model in the infinite 
volume limit is very narrow [l6|. A finite size scaling analysis of the 0(2) spin model's 
magnetization reveals how easy it is to misidentify the universality class of the transition, 
when only small-lattice data is available. 

In section 2 we introduce xQCD. Section 3 introduces the 0(2) spin model and its critical 
behaviour. We describe our Nt = 8 xQCD simulations and present some of our results in 
section 4. In section 5 we present our analysis of the xQCD 'data' using the 0(2) spin 
model. Section 6 is devoted to discussions and conclusions. 
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II. "xQCD" 



In xQCD, the lattice Dirac operator is 

where the first 2 terms are the standard staggered quark Dirac operator for quarks in the 
fundamental representation of colour SU (3), in interaction with the gauge fields on the links 
of the lattice. The auxihary fields a and tt are defined on the sites of the dual lattice and 
the sum is over the 16 sites of the dual lattice closest to the site on which the quark field 
resides. e{x) = The action for these auxihary fields is 

^- = E^^/7(^' + ^') (2) 

where the sum is over the sites of the dual lattice. Since this action is completely local, these 
auxiliary fields can be integrated out introducing a 4-fermion interaction with coupling 
inversely proportional to 7. xQCD preserves the remnant U{1) chiral symmetry of the 
standard staggered quark action for massless quarks. The gauge-field action is the standard 
Wilson action. 

We use hybrid molecular-dynamics simulations to incorporate the factor of [det(A^74)]^J'/^ 
in our simulations. Note that, in contrast to the standard staggered Dirac operator, A^A 
mixes even and odd lattice sites so that we cannot take the square root of the determinant 
by restricting our noisy fermion fields to even (or odd) sites. 

III. THE 0(2) SPIN MODEL IN 3 DIMENSIONS AND CRITICAL BEHAVIOUR 

The simplest version of the 0(2) spin model in 3 dimensions is a non-linear 0(2) sigma 
model. The Hamiltonian H (strictly H/T in condensed matter physics or the action if we 
are studying this model as a field theory) is given by 

H^-JY. Si'Sj-H-Yl^i, (3) 

<i,j> i 

where i and j label sites on a 3 dimensional cubic lattice and < i,j > run over all pairs 
of nearest-neighbour sites. The spins s are 2-vectors of unit length. H is an external, 
symmetry- breaking magnetic field. It is useful to define a temperature T hy T — 1/J. This 
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model and variants have been studied extensively in 



especially useful are listed in references |l2 
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;he literature. Papers we have found 
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The magnetization M of this system is defined as the lattice average of the spins 

1 



M 



V 



(4) 



where V is the volume of the lattice. We first consider the thermodynamic limit V oo. 
At H = 0, in the low temperature phase, 0(2) is broken by a spontaneous magnetization 
with its associated Goldstone mode. At high temperatures the magnetization vanishes and 
0(2) symmetry is restored. The phase transition occurs at a critical point T = = 1/ Jc 
with = 0.454165(4) Q- 

Now let us define those critical exponents and amplitudes which describe the critical 
behaviour of this model and which are relevant to what we might hope to measure. Note 
that these definitions are more general than this specific model. One defines a reduced 
temperature t by t = (T — Tc)/Tc. We first consider the case H = 0. As t ^ 0— the 
magnetization vanishes as 

M = B{-ty^. (5) 



As t 0+, the correlation length 



and the susceptibility 



X = cn- 



Here the susceptibility x is defined by 



X = V{M^). 



(6) 



(7) 



(8) 



The correlation length ^ can be defined in a number of different ways. We will use the 
second moment formula 

" XIF-I 



4sin2(7r/L) 



where L is the length of the box and F is defined by 



Eexp(ip^z)si 



(9) 



(10) 
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where z is the third coordinate of the site i and "Pz = 2tt/L. This form for ^ appears to have 
originated with reference 21[. At small but non-zero H, the magnetization 

M = dH-K (11) 

For the 0(2) spin model these critical exponents have (approximate) values (3m = 0.3490(6), 
z/ = 0.6723(11), 7m ~ 1-32 5 ~ 4.7798, u = 0.79(2) {u is a critical exponent parameterizing 
finite size corrections). 

More information can be obtained about critical points by measuring fluctuation quanti- 
ties (in addition to the susceptibihty) . One such quantity is the fourth-order Binder cumulant 

For if = on an infinite lattice, this quantity is 1, the value for a first-order transition, for 
T < Tc, and 2, the value for a crossover, for T > Tc. Right at Tc it has been measured at 
-B4 = 1.242(2) for the 0(2) spin model [17[. On a finite lattice there is a smooth crossover 
from 1 to 2 as T is varied from small to large values, crossing T^. This crossover steepens as 
V is increased. For large Vs a finite size scaling analysis indicates that these curves cross 
at points which approach [TcB^iTc)] as V ^ 00. Similarly curves of ^/L for finite L (and 
hence V) cross at points which approach Tc and a universal value [^/L]{Tc) = 0.593(2), a 
value specific to 0(2). 

Working at finite volume, where (M) = 0, it is necessary to find an alternative definition 
of the ensemble-averaged magnetization. One which has proved useful is 

M={\M\), (13) 

which will approach the value of the spontaneous magnetization as V ^ 00. It is clear, 
however, that this quantity does not vanish in the high temperature phase except in the 
infinite volume limit. Hence, to test that it has the correct critical behaviour or to measure 
the critical exponent Pm, it is necessary either to extrapolate to infinite volume, or to apply 
a finite size scaling analysis. Finite size scaling for M predicts that close to the critical point 

M(t, L) = L-^ [Q{tL^) + L-'^GitL^) + ...] (14) 

so that plotting M against tL^ should make the magnetization curves for different L 
values coincide, provided L is sufficiently large. 



For QCD, the chiral condensate plays the role of the magnetization. In the 2- dimensional 
space of the remnant U{1) or 0(2) symmetry of staggered fermions, the chiral condensate 
{ipip , i'iplbib'ip) or ((T, tt) corresponds to the magnetization of the 0(2) spin model. 

IV. xQCD SIMULATIONS ON Nt = 8 LATTICES 

We have performed high statistics simulations of massless 2-flavour xQCD on 16'^ x 8 
and 24^ x 8 lattices and lower statistics simulations on 12^ x 8 lattices, using the hybrid 
molecular-dynamics algorithm with dt = 0.05. Measurements are made every 2 1-time-unit 
trajectories. For these simulations the parameter 7 was fixed at 10. At each of 6 values of 
p = 6/^2 piQgg ^he transition (p = 5.5250, 5.5275, 5.5300, 5.5325, 5.5350, 5.5375) we ran 
for 100,000 trajectories on each of the two larger lattices. For (3 = 5.5400, 5.5450, 5.5500 we 
ran for 50,000 trajectories on each of the 2 larger lattices. Shorter runs were made for (3 
values outside of this range. 

Figure Q shows the chiral condensate and Wilson line (Polyakov loop) as functions of 
f3 = Q/ g"^ for both the 16'^ x 8 and the 24'^ x 8 lattices. From this figure, it is clear that the 
finite temperature phase transition occurs somewhere in the range 5.53 (3 5.54. The 
smoothness of the transition suggests strongly that it is second order. The curvature of the 
chiral condensate for the larger lattice (before finite volume rounding sets in) indicates that 
(3m < 1, something that the best finite mass simulations to date have difficulty in observing. 
The reason that such high statistics are needed becomes clearer when one looks at the time 
histories of observables measured during these runs. In figure |2] we show the time history of 
the chiral condensate on the 24^ x 8 lattice for (3 = 5.535 which is close to the transition, 
and (3 = 5.5325 which is just below the transition. Because of critical slowing down, the 
system shows fiuctuations which last for thousands of trajectories close to the critical point. 

With the standard staggered action at zero quark mass, the chiral condensate vanishes 
identically on each configuration, so even if we could simulate at m = we could not 
observe a chiral condensate, even on a configuration by configuration basis. This has to 
be true because the remnant U{1) chiral symmetry is manifest at any fixed value of the 
gauge fields, i.e. on each configuration, which precludes the existence of a chiral condensate. 
With the xQCD action, the condensate is finite on each configuration. This is because, a 
configuration consists of a set of gauge fields and a set of auxiliary fields cr and tt. Any non- 
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FIG. 1: Chiral condensate and Wilson Line as functions of /? = for 16^ x 8 and 24^ x 8 lattices. 



zero set of auxiliary fields breaks chiral symmetry. The vanishing of the chiral condensate 
on a finite volume occurs because, as the system evolves in molecular-dynamics 'time', the 
condensate rotates in the [ipip , iip'j^^^tp) plane so that the ensemble average is zero. This is 
because all U{1) global chiral rotations of the auxiliary fields contribute to the ensemble of 
configurations with equal weight. What we have called (ipip) in the figures above is really 



{Jipi/j"^ ~ V^Ts^s"^^) where the quantities under the square root are lattice averaged before 



being squared. Similarly we define a quantity which we denote (a) which is really {^/o^ + ^l^) . 
In the infinite volume limit each of these quantities will approach the condensate of interest. 
These condensates are related by 

{^Pi;) = 7(a) (15) 



8 



24^x8 Lattice 7=10 jg = 5.535 




20000 40000 60000 80000 100000 
Time = Trajectories 



24^x8 Lattice 7=10 iS = 5.5325 

0.3 I 1 




Time = Trajectories 



FIG. 2: Time histories of the chiral condensate on a 24^ X 8 lattice a) for /3 = 5.535 b) for 
/3 = 5.5325. 
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at dt = 0, (at least in the infinite volume limit) so either quantity can be used as the chiral 
condensate. This will be important when we consider fluctuation quantities, since the ^ipip^ 
we measure is only a stochastic estimator of the condensate, while cr is a true condensate. 
Figureiniis a scatter-plot of the actual chiral condensates measured in the [ipip, ^^'75^5'?/') plane 
on individual conflgurations for 3 (3 values close to the transition, on a X 8 lattice. One 
should flrst note that the condensate does indeed rotate in the plane between conflgurations 
so that the ensemble average is zero, even in the ordered phase. At (3 = 5.5325, below the 
transition, the distribution is concentrated in a ring, indicating a non-zero condensate. By 
(3 = 5.5350, close to the transition, the ring structure is no longer clear, although a peak 
away from the origin still shows up on a radial density proflle. Finally, for f3 = 5.5375, 
above the transition, the maximum density is at the origin indicating that the condensate 
has vanished, and chiral symmetry is restored. 

In the flgure ^ it is clear that to attempt flts to the leading critical scaling behaviour 

(#) = 6(/3,-/5)^'" (3<(3,, (16) 

requires continuation to inflnite volume. Unfortunately, such extrapolations break down 
close to the transition, which is where accurate measurements are most needed. Fits to 
any of these attempted extrapolations always favour values of jSm much smaller than that 
of the 3-dimensional 0(2) spin model. It is for this reason that we have tried fltting to the 
magnetization curves for the 0(2) spin model on small lattices which also incorporate flnite 
volume effects. 

Binder cumulants [2^ have been found useful in characterizing phase transitions. Figured 
shows the 4th order Binder cumulants for the auxiliary flelds {a, n). Following the deflnitions 
given in the previous section, 54 is deflned by 

For low (3, B4 does appear to approach 1 as expected for a flrst order transition (this is 
the flrst order transition when the mass m passes through zero, causing the condensate to 
abruptly change sign). At high (3, i?4 does appear to approach 2, the value for a crossover, as 
expected. In between the curves for the 3 lattice sizes do appear to cross at an intermediate 
value, which is characteristic of a second order transition. The statistics are not good enough 
to determine the positions or values of these intersections, but they appear to be higher than 
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FIG. 3: Scatter-plots of the chiral condensates on a 24^ x 8 lattice. Each of the 50,000 points on 
each graph represents the value on a single configuration the angle brackets here indicate that these 
are configuration averaged quantities but not ensemble averages, a) /3 = 5.5325, b) /3 = 5.5350, c) 
/3 = 5.5375. 

the value 1.242(2) expected if the transition were in the 0(2) universality class. However, 
this probably indicates that our lattices are too small for this intersection to occur very close 
to the infinite volume value. However, this powerful technique should be an excellent way 
of determining the nature of the transitions for larger lattices with very high statistics. 

To get estimates of the position of the transition directly from the 'data' (rather than 
from the fits to spin-model curves) we look at the radial chiral susceptibility, expressed in 
terms of the auxiliary fields 



Xr{cT) = V{{\a\')-{\cT\)'), 

where cr = (a, vr). The position of the phase transition is estimated 
susceptibility, which is found using Ferrenberg-Swendsen reweighting 



(18) 

rom the peak of this 
to interpolate our 
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FIG. 4: The fourth-order binder cumulants for the auxihary fields as functions of (3 for 12^ x 8, 
16^ X 8 and 24^ x 8 lattices. 

'data'. This yields an estimate of Pc = 5.5390(6) for the 16^ x 8 lattice and Pc = 5.5359(5) 
for the 24^ x 8 lattice. (Note that we have also tried using Ferrenberg-Swendsen reweighting 
to interpolate the 'data' for our condensates. While this indicates that the interpolations 
from the measurements at (3 values at the ends of each interval are consistent, the curves it 
yields at current statistics are not sufficiently smooth to help with fits, taking into account 
that the extra points yielded by these interpolations are highly correlated.) In the next 
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section we indicate how (3c calculated in this way approaches the infinite volume limit. 

V. 0(2) SPIN-MODEL ANALYSIS OF xQCD MEASUREMENTS 

We have simulated the 3-dimensional 0(2) spin model with the Hamiltonian given in 
equation 121 on small lattices - 4^, 5'^, 6^, 8^, 9'^, 12'^, 16'^ and 24^. These are adequate for 
the comparisons we make with our results from lattice QCD simulations using the xQCD 
action, which were presented in the previous section. For spin-model results on lattices large 
enough to make comparison with the infinite volume limit, we appeal to the literature. For 
our 4'^ to 12^ lattices, we have interpolated our measurements of the magnetization over the 
range < J < 1.5 using Ferrenberg-Swendsen reweighting. Figure El shows the results of 
our 0(2) simulations and the scaling curves they produced. 

Attempts to perform a finite size scaling analysis (equation I14p around the chiral phase 
transition of xQCD with 0(2) critical indices fail to indicate the existence of a universal 
scaling function. Thus, if the critical point of this theory is in the 0(2) universality class, 
finite volume effects on lattices of the size we use are too large to allow the 'data' to be fit 
by a single scaling function. 

An alternative approach, is to fit the chiral condensates to the magnetization curves for 
the 0(2) spin models on lattices sufficiently small that they too show finite size departures 
from the universal scaling function. Note that these departures depend on the specific form 
of the 0(2) spin- model Hamiltonian used (equation El). We allow ourselves the freedom 
to vary the spin-model lattice sizes to obtain the best fits to the xQCD 'data'. This is 
justified if a single operator is adequate to describe most of the departures from the L — > oo 
universal scaling function of equation El for the 0(2) spin model, and the corresponding 
operator describes most of these departures for xQCD. Varying the lattice size for the 0(2) 
spin model allows us to match the coefficients of the universal scaling function (presumably 
the coefficient G of L^"^ in equation IT^ which describes this departure. If more than one 
non-leading term is required for the model of equation El modeling non-leading behaviour 
of xQCD in this way is somewhat serendipitous. We fit our chiral condensates to the 3- 
parameter form 

i^Pijm = b{M{a{P - Pc) + TJ) (19) 
where {ipip) is the average of the magnitude of the chiral condensate for xQCD and (M) 
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FIG. 5: a) Magnetization as a function of coupling J for the 3-dimensional 
curve is a parametrization of the infinite volume limit given in reference 
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FIG. 6: a) Fit of the chiral condensate for xQCD on a 16^ x 8 lattice to the magnetization of an 

0(2) spin model on a 6^ lattice, b) Fit of the chiral condensate for xQCD on a 24'^ x 8 lattice to 

the magnetization of an 0(2) spin model on an 8'^ lattice. 
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is the magnetization of the 0(2) spin model as a function of T = 1/J. Tc = 1/0.454165 is 
taken as the critical temperature of the 0(2) spin model. jSc the critical value of jS = Q/g"^; a 
and b are the parameters of the fit. For the 16^ x 8 lattice condensate, the best fit is to the 
magnetization on a 6'^ lattice which had a / DOF = 1.36 with parameters jSc = 5.5358(4), 
a = 23.7(8) and b = 0.380(4) over the range 5.50 < /3 < 5.57. Fits to other lattice sizes 
have x^/DOF = 2.34 (4^), x'^/DOF = 1.58 (5^), x^/DOF = 1.55 (8^). For the 24^ x 8 
lattice, the best fits to the condensate are to the magnetizations on lattice sizes 6^ which 
has x^/DOF = 1.85 and parameters = 5.5355(2), a = 53.5(1.3), b = 0.312(2) and 8=^ 
which has x^/DOF = 1.87 with (3^ = 5.5360(2), a = 30.7(6), b = 0.355(2) over the range 
5.51 <P< 5.56. Fits to other lattice sizes have x^lDOF = 2.63 (9^) and x^/DOF = 3.16 
(5^). We are able to obtain reasonable fits over a range of lattice sizes, since if our lattice 
QCD lattices were large enough to be fit by the single scaling function of equation El we 
should be able to obtain good fits of the above form for any 0(2) spin- model lattice large 
enough to also be fit by a single finite-size-scaling function. Figure IHl shows the condensates 
on 16'^ X 8 and 24^ x 8 lattices. The curves are the fits to the 0(2) spin models. We have also 
reanalyzed our 'data' from earlier simulations on an 18'^ x 6 lattice which we had previously 
interpreted as showing tricritical scaling. We find that we can fit these measurements to the 
magnetization for the 0(2) spin model on an 8^ lattice with x^ /DOF = 1.57. Thus A^^^ = 6 
xQCD also shows consistency with 0(2) universality. These fits show the first evidence that 
the chiral condensates are consistent with 0(2) scaling. 

We next examine the correlation length ^, defined as in equation IHl but with cr replacing 
M. One observation that has been made in 0(2) spin models is that curves for C,/L(L is 



the box size) on different (large) lattice sizes cross at a point where ^/L = 0.593(2) 17[, a 
value which is expected to depend only on the universality class. Figure [7| shows C, for the 3 
lattice sizes of our simulations. This gives some indication that the intersection point of the 
correlation lengths is rising towards the universal 0(2) value as the lattice size is increased. 
^/Ns for lattice QCD and ^/L for the spin model should define the same universal scaling 
function in the high temperature phase (which is the only regime where the definition of 
equation El has the interpretation as a correlation length), at corresponding values of (3 and 
T, defined through 

/?=(T-T,)/a + /5„ (20) 
where Pc and a are obtained from the fits of equation Uni except for (a) rather than (iptlj). 
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FIG. 7: The second moment correlation lengths ^ as functions of (3 on 12^ x 8, 16^ x 8 and 24^ x 8 
lattices. Note that = Ng since the spatial 'box' is a cube in each case. 

In figure ISl we have plotted the values of ^/A^s from our xQCD simulations along with ^/L 
from the 0(2) lattice sizes which show the best agreement. The agreement for the larger 
lattice is quite good in the high temperature phase, that for the smaller lattice is only fair. 

The final quantity which we wish to compare with the 0(2) spin model is the chiral 
susceptibility, which is the equivalent of that defined in equation |H1 with M replaced by cr. 
As with the correlation length discussed above, this really only has the interpretation of 
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FIG. 8: (,/N:i for xQCD and ^/L for the 0(2) spin model at corresponding beta values: a) for the 
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susceptibility in the high temperature domain. The values of a and h from fits of (cr) using 
equation Uni give predictions for the chiral susceptibility from the equivalent quantity for the 
spin model 

^^ = V{a^ + 'K^) = Vh'{M'). (21) 

We have plotted these predictions against our xQCD measurements in figure IHl The agree- 
ment between the QCD 'data' and spin-model predictions is very good for the 24'^ x 8 lattice, 
less so for the 16^ x 8 lattice. 

We have used Ferrenberg-Swendsen reweighting to obtain the peaks of the radial suscep- 
tibility for the 0(2) spin model as a function of lattice size, from 4'^ to 32^ lattices. We 
assume that the position of the peaks scales as 

Jpeak = Jc- K L^^ . (22) 

This fit yields = 0.45403(6) from a fit from lattice sizes from 6^ to 32^ compared with 
the best estimate from the literature of Jc = 0.454165(4). The same fit yields 6 = 1.70(2) 
and K = 0.38(2). Hence we conclude that this is a good method for estimating Jc for the 
spin-model and hence should be reasonable for estimating f3c for xQCD. 

We can now study the 0(2) spin model measurements to understand why it is nearly 



V wi 



without 
it was 



impossible to interpret the critical behaviour of our lattice QCD results direct 
comparing them with the spin-model behaviour on small lattices. In reference 
noted that it is only possible to fit the magnetization curve in the infinite volume limit to the 
leading order critical scaling form of equation El over a very small range of \Tc — T\ - roughly 
0-0.07. In terms of J this means from Jc up to perhaps 0.47. Above Jc, the leading finite- 
size corrections to the magnetization are of order 1/L. However, even when J is as large 
as 0.475 extrapolating in 1/L only works for L ^ 12. Hence, for the 5^ to 9^ lattices which 
fit our xQCD 'data', we cannot perform the infinite volume extrapolation reliably. In the 
disordered phase, the leading finite volume corrections are of order 1/VV = 1/L^/^. [This is 
true provided L >> ^, since then one can divide the lattice into N cc V domains of extent 
> ^. The orientations of the spins in these domains will be independent of one another 
so that the magnetization of each spin configuration will be 0{^/N /V) i.e. 0{l/^/V).] 
By the time J is as large as 0.425, departures from this leading order extrapolation are 
already visible on the 9^ lattice. Hence simulations in the region of leading order scaling 
are inaccessible to the small lattices of interest. Even if they were accessible, they represent 
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such a small region of 3 (within ^ 0.0024 of the transition for the 24^ x 8 lattice) as to 

n 

contain only 1 or 2 data points. Reference |16|| finds good fits to the magnetization over a 
much larger range using the scaling form 

M = B{T^ - T)^™[1 + 6i(Tc - T)'^" + h2{T^ - T)] (23) 

- this is the curve plotted in figure El Even if we fixed all 3 critical exponents to their 
0(2) values, this still has 4 parameters. By the time we excluded enough points close to the 
transition and performed the infinite volume extrapolations, such a fit for xQCD would have 
little predictive power. As the authors of this paper and reference j3| note, the narrowness 
of the scaling region appears to be a property peculiar to the 0(2) spin model (at least 
to this particular 0(2) spin model). [Note that the leading departure from scaling (the 
term proportional to hi) in this form corresponds to that in the finite-size scaling form of 
equation IT^ with G{x) — > BbiT^'^x'^'^ as x — > oo. The final term is an order L"^/'^ correction 
term. Although, with the parameters determined in ^| these 2 subleading terms are of 
similar size, it is easy to convince oneself that if the second term is omitted, a reasonable 
approximation to M over the range of interest can be obtained by readjusting parameters 
B and hi] 

Finite size scaling (with just the leading term of equation IT^ apparently fails for the 
chiral condensates measured in our xQCD simulations. It is therefore helpful to apply a 
finite size scaling analysis to the magnetization for 0(2) spin models, including the data 
from small lattices. To do this, we plot the scaling function L^"^^'^M against the scaling 
variable tL^/'^ . Close to the transition {t close to zero), for large enough lattices, the data 
for all lattice sizes should follow the same universal curve. Figure ITOk shows this curve using 
0(2) critical indices (3m and v. Near t = the curves for different lattice sizes coincide. As 
one moves away from t = the curves start to diverge. As the lattice sizes increase the 
curves follow one another for larger ranges of and it appears that they are approaching 

a universal curve. As can be seen, the small lattices do not follow the universal curve very 
far at all. While investigating such scaling, it is interesting to try other values for f3m and 
V. Our earlier, but now suspect, identification of the Nt = 6 transitions as being tricritical, 
suggests we try the critical indices for a tricritical point, jSm = 0.25 and u = 0.5. The 
result is plotted in figure ITUb . The initial impression is that this plot shows much better 
evidence for finite size scaling than that with the correct critical exponents. A closer look 
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reveals that it is the larger lattices that show significant departures from the 'scaling curve'. 
This departure is small but significant for the 12^ lattice, larger for the 16^ lattice and 
embarrassingly large for the 24^ lattice. However, it is clear that if we only had the smaller 
lattices, we would have concluded that the critical point was a tricritical point, not an 0(2) 
critical point. Unfortunately this is essentially the situation for the current round of lattice 
QCD simulations, so we need to be very careful. 



VI. DISCUSSION AND CONCLUSIONS 



We have simulated the finite temperature behaviour of lattice QCD with 2 flavours of 
staggered quarks using the xQCD action which allows us to run at zero quark mass. Run- 
ning with zero quark mass gives us direct access to the second-order phase transition from 
hadronic matter to a quark-gluon plasma, and allows us to to study its universality class. 
Attempts at determining the nature of this transition for staggered quark actions from runs 
at finite mass had led to a state of confusion, since the scaling behaviour near the transition 
often appeared inconsistent with both the 0(4) universality class expected for the contin- 
uum theory and the 0(2) universality class which one would predict from considerations of 
the reduced symmetry of the staggered fermions. The fermion masses that one is forced to 
use with standard staggered actions are too large to clarify what is happening and some 
of the most ambitious simulations of staggered fermion lattice QCD thermodynamics do 
not even try to study the nature of the critical point. Notable exceptions include a recent 
study of this transition at infinite coupling which shows clear evidence for 0(2) universality 



25[, and a study of lattice QCD with 2 light adjoint quark flavours, where the chiral and 



deconfinement transitions are distinct |2^], which also shows evidence for 0(2) universality 
of the chiral transition. 

We present results for the finite temperature behaviour of 2-fiavour lattice QCD with the 
xQCD action on 12'^ x 8, 16^ x 8 and 24^ x 8 lattices. The decrease of the chiral condensate 
as the critical point was approached from below appeared to be too abrupt for 0(2) or 
0(4) universality. However, we found that it was possible to fit the chiral condensates to 
the magnetization curves for an 0(2) spin model on relatively small lattices 5^-9'^. This 
indicates that the transition is consistent with 0(2) universality. The problems with more 
naive approaches to measuring the details of the critical behaviour of this transition is that 
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they do not adequately account for finite volume effects close to the transition, and that for 
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0(2), the non-leading contributions to critical scaling are important except so close to the 
transition as to be inaccessible on the sizes of lattice we use. The spin model predictions 
for the correlation lengths and chiral susceptibilities for lattice QCD are similar enough to 
our measurements to add further support to this 0(2) interpretation. We also made fits to 
the condensates from our 18^ x 6 simulations, which we had earlier interpreted as indicating 
tricritical behaviour, to the 0(2) spin model magnetizations on an 8^ lattice. These fits were 
good enough to convince us that Nt = 6 simulations also yield results consistent with 0(2) 
universality. 

A finite size scaling analysis of our 0(2) spin model magnetizations gave further clues as 
to why the behaviour of our lattice QCD measurements hinted at a tricritical explanation. 
The approach of the 0(2) spin model magnetizations to the universal finite size scaling curve 
with increasing lattice size is slow. However, if we perform a finite size scaling analysis using 
the critical exponents of the tricritical point, rather than those appropriate for 0(2), we find 
a rapid approach to an apparently universal curve with increasing lattice size. It is only 
by going to larger lattices, where the measurements start to diverge from this false scaling 
curve, that the tricritical interpretation is invalidated. For lattice QCD we do not at present 
have the luxury of being able to perform simulations on much larger lattices. 

The results we presented in this paper are for a single small fixed value (in lattice units) 
of the 4-fermion coupling. It would be helpful if we could perform simulations at another 
(even smaller) value of this coupling, since the critical value of /3, 5.535 ^ /5c ~ 5.536 is 
significantly above that for the standard staggered action at Nt = 8, where reference j^l 
estimates Pc = 5.44(3). We should also extend our analyses to finite quark masses to study 
critical scaling involving the critical exponent S. We note that we could extend our current 
measurements to finite quark masses using Ferrenberg-Swendsen extrapolation. However, 
our attempts to do such an extrapolation indicated that this is only reliable for a range of 
masses where the condensate vanishes linearly with decreasing quark mass. For finite mass 
studies one needs quark masses large enough to get beyond this linear regime. 

So far, we have only compared the critical behaviour of lattice QCD with the 0(2) spin 
model. We need to extend our studies to spin models in other universality classes as for 
example the 0(4) or the mean- field universality classes, before we can say with any certainty 
whether the evidence for 0(2) universality is compelling. Such studies are planned for the 
near future, since the spin-model simulations they involve are relatively inexpensive. 
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Finally, we need to address finite dt effects. We have found that finite dt effects can 
change the apparent nature of the phase transition in our studies of QCD at finite isospin 



chemical potential [28|, |29| (see also 



30|). There the problem was that the difference between 



the effective /3, calculated by the equipartition theorem from the molecular-dynamics kinetic 
energies and the input /?, a finite dt effect, was much larger below the transition than above, 
which distorted the transition. We have measured such (3 shifts in our xQCD simulations, 
and while they are as large as those measured in our QCD at finite isospin simulations, they 
do not significantly change as we cross the transition, and are therefore unlikely to distort 
this transition. Nevertheless we intend to convert our simulations from hybrid molecular- 

n 

dynamics to exact RHMC simulations [31| . Because the lower bound of the spectrum of the 
xQCD Dirac operator is unknown, this is not an entirely trivial conversion. However, our 
experience with using the RHMC method to simulate QCD at finite isospin density, where 
the lower bound on the Dira c op erator is also unknown, indicates that RHMC simulations 
can be applied in such cases j^^l . 
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